DEseq2 Analysis

Install required packages

if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install("DESeq2")
install.packages("pheatmap")
install.packages("RColorBrewer")
install.packages("dplyr")
install.packages("ggplot2")
BiocManager::install("genefilter")
BiocManager::install("ashr")
BiocManager::install("VennDetail")
install.packages("tidyverse")
install.packages("kableExtra")

Load required packages

library(DESeq2)
library(RColorBrewer)
library(dplyr)
library(ggplot2)
library(genefilter)
library(ashr)
library(VennDetail)
library(tidyverse)
library(kableExtra)
library(pheatmap)
library(ComplexHeatmap)
library(circlize)

Count Matrix Construction

pheno.data <- read.csv("pheno_data.csv")
gene.count.matrix <- read.table(file = "gene_count_matrix.csv",header = T,
                                sep = ",",row.names = NULL)
head(gene.count.matrix)

##                   gene_id sample_01 sample_02 sample_03 sample_04 sample_05
## 1      Cre14.g619950.v5.5       115       127        93       162       167
## 2      Cre01.g048900.v5.5         0         0        12        10        42
## 3 Cre12.g543400.v5.5|FDH1       131       110       114       118       175
## 4      Cre16.g685837.v5.5        63        56        53        29       109
## 5      Cre02.g119526.v5.5         0         0         0         0         0
## 6      Cre06.g278159.v5.5      2707      2417      2316      2354      3916
##   sample_06 sample_07 sample_08 sample_09 sample_10 sample_11 sample_12
## 1       127       128       188       134        66        81        88
## 2        36        28         9        33        17        28        22
## 3       138       180       170        62        92        70        76
## 4       108        89       115        55        68        51        45
## 5         0         0         0         0         0         0         0
## 6      3663      3556      3537      1918      1959      1914      1826
##   sample_13 sample_14 sample_15 sample_16 sample_17 sample_18 sample_19
## 1        21         9        24        22        67        73        59
## 2         0        11        10         0        60        47        56
## 3        44        12        47        34       197       169       175
## 4        27        35        20        19        94       117       104
## 5         0         0         0         0         0         0         0
## 6       737       807       649       688      3347      3055      2937
##   sample_20 sample_21 sample_22 sample_23 sample_24 sample_25 sample_26
## 1        69        76        73        52        51       143        99
## 2        39        22        15        24         0        57        50
## 3       167       175       233       223       211        52        66
## 4        89        45        56        50        36        83        65
## 5         0         0         0         0         0         0         0
## 6      3173      2702      2480      2466      2429      2193      2079
##   sample_27 sample_28 sample_29 sample_30 sample_31 sample_32 sample_33
## 1       136        93       143       122       122       158       166
## 2        54        43         0         0         0         0        28
## 3        77        58        96       131        98       121       170
## 4        56        47       107        78       110        78       107
## 5         0         0         0         0         0         0         0
## 6      2104      1966      4211      4262      4106      3964      3586
##   sample_34 sample_35 sample_36 sample_37 sample_38 sample_39 sample_40
## 1       175       145       105        48        46        44        49
## 2        17        11        21         0         7         0         0
## 3       146       186        81        87        61        63        71
## 4        97       114       103        30        65        59        31
## 5         0         0         0         0         0         0         0
## 6      3233      3417      3375      1573      1519      1486      1466
##   sample_41 sample_42 sample_43 sample_44 sample_45 sample_46 sample_47
## 1       111        81        61        63         0         0         0
## 2        72        63        31        77         0         0         0
## 3       205       167       158       151        91       106        82
## 4       109        61        61        88        62        61        63
## 5         0         0         0         0         0         0         0
## 6      4477      4189      4127      4133      2854      2824      2724
##   sample_48
## 1         5
## 2         0
## 3       101
## 4        58
## 5         0
## 6      2738

gene.ids <- sapply(X = strsplit(x = gene.count.matrix$gene_id,split = "\\|"),
                   FUN = function(x){return(x[1])})
gene.count.matrix <- gene.count.matrix[,-1]

rownames(gene.count.matrix) = make.names(gene.ids, unique=TRUE)
head(gene.count.matrix)

##                    sample_01 sample_02 sample_03 sample_04 sample_05 sample_06
## Cre14.g619950.v5.5       115       127        93       162       167       127
## Cre01.g048900.v5.5         0         0        12        10        42        36
## Cre12.g543400.v5.5       131       110       114       118       175       138
## Cre16.g685837.v5.5        63        56        53        29       109       108
## Cre02.g119526.v5.5         0         0         0         0         0         0
## Cre06.g278159.v5.5      2707      2417      2316      2354      3916      3663
##                    sample_07 sample_08 sample_09 sample_10 sample_11 sample_12
## Cre14.g619950.v5.5       128       188       134        66        81        88
## Cre01.g048900.v5.5        28         9        33        17        28        22
## Cre12.g543400.v5.5       180       170        62        92        70        76
## Cre16.g685837.v5.5        89       115        55        68        51        45
## Cre02.g119526.v5.5         0         0         0         0         0         0
## Cre06.g278159.v5.5      3556      3537      1918      1959      1914      1826
##                    sample_13 sample_14 sample_15 sample_16 sample_17 sample_18
## Cre14.g619950.v5.5        21         9        24        22        67        73
## Cre01.g048900.v5.5         0        11        10         0        60        47
## Cre12.g543400.v5.5        44        12        47        34       197       169
## Cre16.g685837.v5.5        27        35        20        19        94       117
## Cre02.g119526.v5.5         0         0         0         0         0         0
## Cre06.g278159.v5.5       737       807       649       688      3347      3055
##                    sample_19 sample_20 sample_21 sample_22 sample_23 sample_24
## Cre14.g619950.v5.5        59        69        76        73        52        51
## Cre01.g048900.v5.5        56        39        22        15        24         0
## Cre12.g543400.v5.5       175       167       175       233       223       211
## Cre16.g685837.v5.5       104        89        45        56        50        36
## Cre02.g119526.v5.5         0         0         0         0         0         0
## Cre06.g278159.v5.5      2937      3173      2702      2480      2466      2429
##                    sample_25 sample_26 sample_27 sample_28 sample_29 sample_30
## Cre14.g619950.v5.5       143        99       136        93       143       122
## Cre01.g048900.v5.5        57        50        54        43         0         0
## Cre12.g543400.v5.5        52        66        77        58        96       131
## Cre16.g685837.v5.5        83        65        56        47       107        78
## Cre02.g119526.v5.5         0         0         0         0         0         0
## Cre06.g278159.v5.5      2193      2079      2104      1966      4211      4262
##                    sample_31 sample_32 sample_33 sample_34 sample_35 sample_36
## Cre14.g619950.v5.5       122       158       166       175       145       105
## Cre01.g048900.v5.5         0         0        28        17        11        21
## Cre12.g543400.v5.5        98       121       170       146       186        81
## Cre16.g685837.v5.5       110        78       107        97       114       103
## Cre02.g119526.v5.5         0         0         0         0         0         0
## Cre06.g278159.v5.5      4106      3964      3586      3233      3417      3375
##                    sample_37 sample_38 sample_39 sample_40 sample_41 sample_42
## Cre14.g619950.v5.5        48        46        44        49       111        81
## Cre01.g048900.v5.5         0         7         0         0        72        63
## Cre12.g543400.v5.5        87        61        63        71       205       167
## Cre16.g685837.v5.5        30        65        59        31       109        61
## Cre02.g119526.v5.5         0         0         0         0         0         0
## Cre06.g278159.v5.5      1573      1519      1486      1466      4477      4189
##                    sample_43 sample_44 sample_45 sample_46 sample_47 sample_48
## Cre14.g619950.v5.5        61        63         0         0         0         5
## Cre01.g048900.v5.5        31        77         0         0         0         0
## Cre12.g543400.v5.5       158       151        91       106        82       101
## Cre16.g685837.v5.5        61        88        62        61        63        58
## Cre02.g119526.v5.5         0         0         0         0         0         0
## Cre06.g278159.v5.5      4127      4133      2854      2824      2724      2738

Creating DEseq data set object

dds <- DESeqDataSetFromMatrix(countData=gene.count.matrix, colData=pheno.data, 
                              design = ~genotype+ condition +  
                                genotype:condition)

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors

dds$genotype

##  [1] wt   wt   wt   wt   wt   wt   wt   wt   wt   wt   wt   wt   atg8 atg8 atg8
## [16] atg8 atg8 atg8 atg8 atg8 atg8 atg8 atg8 atg8 wt   wt   wt   wt   wt   wt  
## [31] wt   wt   wt   wt   wt   wt   atg8 atg8 atg8 atg8 atg8 atg8 atg8 atg8 atg8
## [46] atg8 atg8 atg8
## Levels: atg8 wt

dds$condition

##  [1] 0h 0h 0h 0h 0h 0h 0h 0h 0h 0h 0h 0h 0h 0h 0h 0h 0h 0h 0h 0h 0h 0h 0h 0h 4h
## [26] 4h 4h 4h 4h 4h 4h 4h 4h 4h 4h 4h 4h 4h 4h 4h 4h 4h 4h 4h 4h 4h 4h 4h
## Levels: 0h 4h

dds

## class: DESeqDataSet 
## dim: 17749 48 
## metadata(1): version
## assays(1): counts
## rownames(17749): Cre14.g619950.v5.5 Cre01.g048900.v5.5 ...
##   Cre21.g753197.v5.5 Cre10.g441750.v5.5
## rowData names(0):
## colnames(48): sample_01 sample_02 ... sample_47 sample_48
## colData names(3): sample genotype condition

nrow(dds)

## [1] 17749

# Our experiment consists of 12 samples each with 3 replicates:

dds$sample <- factor(paste0("sample_", rep(c(1,2,3,4,5,6,7,8,9,10,11,12), 
                                           each=4))
                     , levels = c("sample_1", "sample_2","sample_3","sample_4",
                                  "sample_5","sample_6",
                                  "sample_7","sample_8","sample_9","sample_10",
                                  "sample_11","sample_12"))
dds$sample

##  [1] sample_1  sample_1  sample_1  sample_1  sample_2  sample_2  sample_2 
##  [8] sample_2  sample_3  sample_3  sample_3  sample_3  sample_4  sample_4 
## [15] sample_4  sample_4  sample_5  sample_5  sample_5  sample_5  sample_6 
## [22] sample_6  sample_6  sample_6  sample_7  sample_7  sample_7  sample_7 
## [29] sample_8  sample_8  sample_8  sample_8  sample_9  sample_9  sample_9 
## [36] sample_9  sample_10 sample_10 sample_10 sample_10 sample_11 sample_11
## [43] sample_11 sample_11 sample_12 sample_12 sample_12 sample_12
## 12 Levels: sample_1 sample_2 sample_3 sample_4 sample_5 sample_6 ... sample_12

dds <- collapseReplicates(dds, dds$sample)

# Only keep counts without 0:

keep <- rowSums(counts(dds) > 1) >=3
dds <- dds[keep,]
nrow(dds)

## [1] 13977

gene.ids <- rownames(dds)

Set WT and condition 0h as the reference level

dds$genotype <- relevel(dds$genotype, ref = "wt")
dds$condition <- relevel(dds$condition, ref = "0h")

Executing DEseq2

dds <- DESeq(dds, minReplicatesForReplace = Inf) 
resultsNames(dds) 

## [1] "Intercept"                "genotype_atg8_vs_wt"     
## [3] "condition_4h_vs_0h"       "genotypeatg8.condition4h"

# Data visualization according to Cooks coefficient

par(mar=c(8,5,2,2))
boxplot(log10(assays(dds)[["cooks"]]), range=0, las=2, col=rainbow(16)) 

Data previsualization

Data transformations

Normal logarithm, variance stabilizing and regularized transformations were performed in order to visualize data; finally, the variance stabilizing transformation was selected.

ntd <- normTransform(dds)
vsd <- vst(dds, blind=FALSE)

## -- note: fitType='parametric', but the dispersion trend was not well captured by the
##    function: y = a/x + b, and a local regression fit was automatically substituted.
##    specify fitType='local' or 'mean' to avoid this message next time.

rld <- rlog(dds, blind = FALSE)

df <- bind_rows(as_data_frame(assay(ntd)[, 1:2]) %>% mutate(transformation = "ntd"),
                as_data_frame(assay(vsd)[, 1:2]) %>% mutate(transformation = "vst"),
                as_data_frame(assay(rld)[, 1:2]) %>% mutate(transformation = "rlog"))

## Warning: `as_data_frame()` was deprecated in tibble 2.0.0.
## Please use `as_tibble()` instead.
## The signature and semantics have changed, see `?as_tibble`.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.

# Plot exploring every transformation

colnames(df)[1:2] <- c("x", "y")  
lvls <- c("ntd", "vst", "rlog")
df$transformation <- factor(df$transformation, levels=lvls)
ggplot(df, aes(x = x, y = y)) + geom_hex(bins = 80) +
  coord_fixed() + facet_grid( . ~ transformation) 

Principal Component Analysis

pcaData <- plotPCA(vsd, intgroup=c("genotype", "condition"), returnData=TRUE)
percentVar <- round(100 * attr(pcaData, "percentVar"))
ggplot(pcaData, aes(PC1, PC2, color=genotype, shape=condition)) +
  geom_point(size=3) +
  xlab(paste0("PC1: ",percentVar[1],"% variance")) +
  ylab(paste0("PC2: ",percentVar[2],"% variance")) + 
  coord_fixed()

Gene clustering visualization

topVarGenes <- head(order(rowVars(assay(vsd)), decreasing = TRUE), 20)
mat  <- assay(vsd)[ topVarGenes, ]
mat  <- mat - rowMeans(mat)
anno <- as.data.frame(colData(vsd)[, c("genotype","condition")])
pheatmap(mat, annotation_col = anno)

Gene Differential Expression Analysis

Five different contrasts were performed in order to understand the effects of cerulenin in C. reinhardtii cells. Those were:

• Effect of cerulenin in wild-type (0h vs 4h)
• Effect of cerulenin in mutant (0h vs 4h)
• Difference between wild-type and mutant 0h
• Difference between wild-type and mutant 4h
• Different response for both genotypes (the interaction term)

Load some useful markers

ATG_markers <- read.table(file = "ATG marker genes .txt",header = T,
                                sep = "\t",row.names = NULL)

Chloro_stress_markers <- read.table(file = "Chloroplast chaperones and stress marker genes.txt",header =T, 
                                    sep = "\t",row.names = NULL)

Chloro_redox_markers <- read.table(file = "Choloplast redox marker genes.txt",header = T,
                                sep = "\t",row.names = NULL)

Effect of cerulenin in wild-type (0h vs 4h)

Gene obtention

resWT = results(dds, contrast=c("condition","4h","0h"), cooksCutoff = FALSE, 
                independentFiltering = FALSE)
resWT

## log2 fold change (MLE): condition 4h vs 0h 
## Wald test p-value: condition 4h vs 0h 
## DataFrame with 13977 rows and 6 columns
##                      baseMean log2FoldChange     lfcSE       stat    pvalue
##                     <numeric>      <numeric> <numeric>  <numeric> <numeric>
## Cre14.g619950.v5.5   307.8280      0.0281565  0.717068  0.0392662  0.968678
## Cre01.g048900.v5.5    76.4131      0.3799222  1.828657  0.2077602  0.835416
## Cre12.g543400.v5.5   428.9410     -0.2585408  0.320513 -0.8066479  0.419869
## Cre16.g685837.v5.5   255.7288      0.2418724  0.453962  0.5328032  0.594170
## Cre06.g278159.v5.5  9770.0557      0.1504420  0.124149  1.2117907  0.225592
## ...                       ...            ...       ...        ...       ...
## Cre06.g278104.v5.5  889.54623     -0.4304988  0.309745 -1.3898503 0.1645743
## Cre17.g717150.v5.5  146.81853      0.0579478  1.170897  0.0494901 0.9605287
## Cre10.g443600.v5.5  170.10265     -1.9623294  1.081889 -1.8137987 0.0697087
## Cre14.g621400.v5.5    5.13639     -0.1021788  4.276178 -0.0238949 0.9809364
## Cre10.g441750.v5.5 1599.45474      0.1400955  0.181245  0.7729615 0.4395452
##                         padj
##                    <numeric>
## Cre14.g619950.v5.5  0.997405
## Cre01.g048900.v5.5  0.990087
## Cre12.g543400.v5.5  0.862068
## Cre16.g685837.v5.5  0.942215
## Cre06.g278159.v5.5  0.652547
## ...                      ...
## Cre06.g278104.v5.5  0.550696
## Cre17.g717150.v5.5  0.996460
## Cre10.g443600.v5.5  0.324773
## Cre14.g621400.v5.5  0.997405
## Cre10.g441750.v5.5  0.870623

summary(resWT)

## 
## out of 13977 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0 (up)       : 909, 6.5%
## LFC < 0 (down)     : 980, 7%
## outliers [1]       : 0, 0%
## low counts [2]     : 0, 0%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

log.fold.change_WT <- resWT$log2FoldChange
q.value_WT <- resWT$padj
names(log.fold.change_WT) <- gene.ids
names(q.value_WT) <- gene.ids
base.meanWT <- resWT$baseMean
names(base.meanWT) <- gene.ids

# Take a look at upregulated and downregulated genes: 

top_20_overexpressed_genes_WT <- resWT[order(resWT$log2FoldChange, decreasing = T),]
kable(top_20_overexpressed_genes_WT[1:20,-(3:4)])

	baseMean	log2FoldChange	pvalue	padj
Cre10.g465763.v5.5	17.445204	20.612283	0.0000012	0.0000297
Cre07.g331350.v5.5	19.178658	20.226262	0.0000004	0.0000097
Cre10.g466200.v5.5	15.235117	18.706363	0.0000109	0.0002224
Cre02.g143107.v5.5	7.114070	17.784866	0.0000295	0.0005527
Cre06.g250650.v5.5	13.336058	15.503286	0.0002575	0.0037533
Cre02.g102300.v5.5	33.165537	8.911089	0.0019320	0.0205662
Cre17.g735876.v5.5	22.382807	8.046166	0.0126038	0.0945220
Cre02.g116900.v5.5	12.062879	7.629730	0.0605932	0.2956062
Cre08.g364850.v5.5	8.610772	7.489975	0.0320338	0.1896386
Cre16.g680454.v5.5	350.851311	7.410076	0.0001349	0.0021043
Cre02.g091226.v5.5	16.082502	7.296935	0.0859688	0.3740927
Cre12.g511952.v5.5	11.269179	7.184714	0.0909203	0.3867668
Cre02.g095151.v5.5	13189.468430	7.063199	0.0000000	0.0000000
Cre07.g341750.v5.5	36.286173	6.880834	0.0077292	0.0644578
Cre08.g369200.v5.5	5.101397	6.578353	0.1158762	0.4496394
Cre12.g540351.v5.5	46.825918	6.539815	0.0033205	0.0323643
Cre03.g167690.v5.5	3.495144	6.501141	0.0738629	0.3369917
Cre03.g144444.v5.5	8.399966	6.444530	0.1295947	0.4817355
Cre04.g229422.v5.5	3.341127	6.298153	0.0952044	0.3978966
Cre12.g554929.v5.5	123.196764	6.295387	0.0000011	0.0000260

top_20_overrepressed_genes_WT <- resWT[order(resWT$log2FoldChange, decreasing = F),]
kable(top_20_overrepressed_genes_WT[1:20,-(3:4)])

	baseMean	log2FoldChange	pvalue	padj
Cre01.g043250.v5.5	102.22785	-22.99906	0e+00	0.0e+00
Cre16.g658300.v5.5	46.21710	-22.89138	0e+00	0.0e+00
Cre02.g115400.v5.5	20.74418	-22.35411	0e+00	9.0e-07
Cre12.g500250.v5.5	28.14109	-22.03587	0e+00	0.0e+00
Cre03.g179400.v5.5	39.29917	-22.00429	0e+00	0.0e+00
Cre06.g292249.v5.5	18.96145	-21.81929	1e-07	2.8e-06
Cre06.g283000.v5.5	20.54063	-21.80617	0e+00	1.0e-07
Cre17.g735750.v5.5	27.85918	-21.60836	0e+00	0.0e+00
Cre01.g034380.v5.5	38.57187	-21.42391	0e+00	0.0e+00
Cre13.g588250.v5.5	49.57133	-21.38082	0e+00	0.0e+00
Cre16.g684715.v5.5	14.85207	-21.37278	4e-07	9.9e-06
Cre03.g143827.v5.5	66.50156	-21.31874	0e+00	0.0e+00
Cre14.g629840.v5.5	37.01919	-21.31223	0e+00	0.0e+00
Cre14.g632350.v5.5	93.48560	-21.17466	0e+00	0.0e+00
Cre12.g554100.v5.5	20.57133	-21.16066	0e+00	2.0e-07
Cre09.g402050.v5.5	54.88133	-21.12392	0e+00	0.0e+00
Cre02.g083050.v5.5	18.89988	-21.07576	1e-07	4.0e-06
Cre09.g397050.v5.5	47.73692	-20.98995	0e+00	0.0e+00
Cre16.g658400.v5.5	59.75973	-20.89439	0e+00	0.0e+00
Cre04.g226550.v5.5	12.19448	-20.86036	9e-07	2.3e-05

# Differentialy activated and repressed genes with a log fold change and p-adj threshold of 1.5 and 0.05

activated.genes.deseq2_WT <- gene.ids[log.fold.change_WT > 1 & q.value_WT < 0.05]
activated.genes.deseq2_WT <- activated.genes.deseq2_WT[!is.na(activated.genes.deseq2_WT)]
activated.genes.deseq2_WT_good <- substr(activated.genes.deseq2_WT,1,nchar(activated.genes.deseq2_WT)-5)
write.table(activated.genes.deseq2_WT_good,file="activated.genes.deseq2_WT.txt",sep="\t",quote=F,row.names = F)

repressed.genes.deseq2_WT <- gene.ids[log.fold.change_WT < - 1 & q.value_WT < 0.05]
repressed.genes.deseq2_WT <- repressed.genes.deseq2_WT[!is.na(repressed.genes.deseq2_WT)]
repressed.genes.deseq2_WT_good <- substr(repressed.genes.deseq2_WT,1,nchar(repressed.genes.deseq2_WT)-5)
write.table(repressed.genes.deseq2_WT_good,file="repressed.genes.deseq2_WT.txt",sep="\t",quote=F,row.names = F)

length(activated.genes.deseq2_WT_good)

## [1] 458

length(repressed.genes.deseq2_WT_good)

## [1] 350

MA plot

plotMA(resWT,alpha=0.05, ylim=c(-30,30), colSig="grey")
abline(h=c(-1,1), col="black", lwd=1)
points(x = base.meanWT[activated.genes.deseq2_WT],
       y = log.fold.change_WT[activated.genes.deseq2_WT],col="red",cex=0.2,pch=19)
points(x = base.meanWT[repressed.genes.deseq2_WT],
       y = log.fold.change_WT[repressed.genes.deseq2_WT],col="blue",cex=0.2,pch=19)
text(x =100000, y = 15, label = length(activated.genes.deseq2_WT_good),
     col = "red",   # Color del texto
     font = 2,      # Estilo
     cex = 1.5)     # Tamaño
text(x =100000, y = -15, label = length(repressed.genes.deseq2_WT_good),
     col = "blue",   # Color del texto
     font = 2,      # Estilo
     cex = 1.5)     # Tamaño

legend("bottomright", legend = c("Overexpressed", "Repressed"),
       lwd = 3, col = c("red", "blue"))

Effect of cerulenin in mutant (0h vs 4h)

Gene obtention

res_atg8 <- results(dds, list(c("genotypeatg8.condition4h","condition_4h_vs_0h")),
                    cooksCutoff = FALSE, independentFiltering = FALSE)
res_atg8

## log2 fold change (MLE): genotypeatg8.condition4h+condition_4h_vs_0h effect 
## Wald test p-value: genotypeatg8.condition4h+condition_4h_vs_0h effect 
## DataFrame with 13977 rows and 6 columns
##                      baseMean log2FoldChange     lfcSE       stat      pvalue
##                     <numeric>      <numeric> <numeric>  <numeric>   <numeric>
## Cre14.g619950.v5.5   307.8280     -0.0177850  0.721654 -0.0246447 9.80338e-01
## Cre01.g048900.v5.5    76.4131     -0.3532859  1.831528 -0.1928914 8.47044e-01
## Cre12.g543400.v5.5   428.9410     -0.0527444  0.321006 -0.1643096 8.69487e-01
## Cre16.g685837.v5.5   255.7288      0.3030096  0.456488  0.6637846 5.06828e-01
## Cre06.g278159.v5.5  9770.0557      0.5915680  0.124397  4.7554709 1.97984e-06
## ...                       ...            ...       ...        ...         ...
## Cre06.g278104.v5.5  889.54623       0.402382  0.311428   1.292054    0.196338
## Cre17.g717150.v5.5  146.81853      -1.580214  1.173164  -1.346968    0.177991
## Cre10.g443600.v5.5  170.10265      -0.202010  1.079760  -0.187088    0.851592
## Cre14.g621400.v5.5    5.13639      -2.283619  3.977653  -0.574112    0.565892
## Cre10.g441750.v5.5 1599.45474       0.148161  0.181933   0.814375    0.415430
##                           padj
##                      <numeric>
## Cre14.g619950.v5.5 9.98339e-01
## Cre01.g048900.v5.5 9.70421e-01
## Cre12.g543400.v5.5 9.75746e-01
## Cre16.g685837.v5.5 8.13498e-01
## Cre06.g278159.v5.5 3.98446e-05
## ...                        ...
## Cre06.g278104.v5.5    0.517192
## Cre17.g717150.v5.5    0.490686
## Cre10.g443600.v5.5    0.972125
## Cre14.g621400.v5.5    0.847746
## Cre10.g441750.v5.5    0.747101

log.fold.change_atg8 <- res_atg8$log2FoldChange
q.value_atg8 <- res_atg8$padj
names(log.fold.change_atg8) <- gene.ids
names(q.value_atg8) <- gene.ids
base.mean_atg8 <- res_atg8$baseMean
names(base.mean_atg8) <- gene.ids

# Take a look at upregulated and downregulated genes: 

top_20_overexpressed_genes_atg8 <- res_atg8[order(res_atg8$log2FoldChange, decreasing = T),]
kable(top_20_overexpressed_genes_atg8[1:20,-(3:4)])

	baseMean	log2FoldChange	pvalue	padj
Cre01.g031700.v5.5	71.40996	22.31599	0e+00	0.00e+00
Cre03.g169150.v5.5	84.33879	22.21561	0e+00	0.00e+00
Cre03.g206900.v5.5	39.26690	22.20814	0e+00	0.00e+00
Cre05.g242750.v5.5	31.80983	22.17803	0e+00	0.00e+00
Cre03.g183150.v5.5	52.34164	22.08262	0e+00	0.00e+00
Cre16.g650000.v5.5	31.64993	21.81398	0e+00	2.00e-07
Cre01.g055457.v5.5	18.41569	21.70069	3e-07	8.00e-06
Cre06.g282826.v5.5	42.54245	21.65897	0e+00	0.00e+00
Cre16.g688900.v5.5	26.82071	21.59666	0e+00	3.00e-07
Cre24.g755097.v5.5	25.65882	21.46373	0e+00	4.00e-07
Cre06.g305302.v5.5	18.84289	21.35537	5e-07	1.18e-05
Cre07.g346900.v5.5	87.77094	21.16046	0e+00	0.00e+00
Cre01.g009650.v5.5	19.95885	21.02765	0e+00	2.00e-07
Cre13.g581700.v5.5	50.68467	20.99081	0e+00	0.00e+00
Cre13.g569000.v5.5	31.53099	20.90911	0e+00	0.00e+00
Cre15.g637315.v5.5	10.27450	20.87907	9e-07	1.98e-05
Cre12.g544113.v5.5	21.59394	20.87294	0e+00	2.00e-07
Cre17.g724873.v5.5	12.39615	20.76123	1e-06	2.25e-05
Cre04.g216700.v5.5	20.06815	20.75520	0e+00	3.00e-07
Cre02.g109900.v5.5	70.00800	20.72860	0e+00	0.00e+00

top_20_overrepressed_genes_atg8 <- res_atg8[order(res_atg8$log2FoldChange, decreasing = F),]
kable(top_20_overrepressed_genes_atg8[1:20,-(3:4)])

	baseMean	log2FoldChange	pvalue	padj
Cre06.g283000.v5.5	20.54063	-29.73628	0e+00	0.0e+00
Cre17.g735750.v5.5	27.85918	-29.05778	0e+00	0.0e+00
Cre06.g292249.v5.5	18.96145	-26.98877	0e+00	0.0e+00
Cre12.g500250.v5.5	28.14109	-26.72339	0e+00	0.0e+00
Cre02.g115400.v5.5	20.74418	-25.68656	0e+00	0.0e+00
Cre16.g684715.v5.5	14.85207	-25.34339	0e+00	1.0e-07
Cre06.g294250.v5.5	123.06077	-24.04041	0e+00	0.0e+00
Cre08.g376740.v5.5	77.07946	-23.39426	0e+00	0.0e+00
Cre03.g205800.v5.5	69.35132	-23.21297	0e+00	1.3e-06
Cre13.g604550.v5.5	56.05740	-23.16705	0e+00	0.0e+00
Cre07.g315050.v5.5	88.92973	-23.13839	0e+00	0.0e+00
Cre09.g412450.v5.5	61.27340	-23.03721	0e+00	0.0e+00
Cre10.g428000.v5.5	52.12164	-23.01926	0e+00	0.0e+00
Cre13.g592350.v5.5	49.26531	-22.72924	0e+00	0.0e+00
Cre07.g325713.v5.5	47.76235	-22.60765	0e+00	0.0e+00
Cre06.g278168.v5.5	57.62603	-22.41235	0e+00	0.0e+00
Cre12.g547727.v5.5	29.91381	-22.32585	0e+00	0.0e+00
Cre05.g246250.v5.5	43.36180	-22.26541	0e+00	0.0e+00
Cre01.g010832.v5.5	37.26769	-22.11360	0e+00	0.0e+00
Cre16.g695200.v5.5	28.69367	-22.10441	1e-07	2.5e-06

activated.genes.deseq2_atg8 <- gene.ids[log.fold.change_atg8 > 1 & q.value_atg8 < 0.05]
activated.genes.deseq2_atg8 <- activated.genes.deseq2_atg8[!is.na(activated.genes.deseq2_atg8)]
activated.genes.deseq2_atg8_good <- substr(activated.genes.deseq2_atg8,1,nchar(activated.genes.deseq2_atg8)-5)
write.table(activated.genes.deseq2_atg8_good,file="activated.genes.deseq2_atg8.txt",sep="\t",quote=F,row.names = F)

repressed.genes.deseq2_atg8 <- gene.ids[log.fold.change_atg8 < -1 & q.value_atg8 < 0.05]
repressed.genes.deseq2_atg8 <- repressed.genes.deseq2_atg8[!is.na(repressed.genes.deseq2_atg8)]
repressed.genes.deseq2_atg8_good <- substr(repressed.genes.deseq2_atg8,1,nchar(repressed.genes.deseq2_atg8)-5)
write.table(repressed.genes.deseq2_atg8_good,file="repressed.genes.deseq2_atg8.txt",sep="\t",quote=F,row.names = F)

length(activated.genes.deseq2_atg8_good)

## [1] 618

length(repressed.genes.deseq2_atg8_good)

## [1] 441

MA plot

plotMA(res_atg8,alpha=0.05, ylim=c(-30,30),colSig="grey")
abline(h=c(-1,1), col="black", lwd=1)
points(x = base.meanWT[activated.genes.deseq2_atg8],
       y = log.fold.change_atg8[activated.genes.deseq2_atg8],col="red",cex=0.2,pch=19)
points(x = base.mean_atg8[repressed.genes.deseq2_atg8],
       y = log.fold.change_atg8[repressed.genes.deseq2_atg8],col="blue",cex=0.2,pch=19)
text(x =100000, y = 15, label = length(activated.genes.deseq2_atg8_good),
     col = "red",   # Color del texto
     font = 2,      # Estilo
     cex = 1.5)     # Tamaño
text(x =100000, y = -15, label = length(repressed.genes.deseq2_atg8_good),
     col = "blue",   # Color del texto
     font = 2,      # Estilo
     cex = 1.5)     # Tamaño

legend("bottomright", legend = c("Overexpressed", "Repressed"),
       lwd = 3, col = c("red", "blue"))

Plotting log2 fold change with target genes for wt and atg8 treated with cerulenin

# Test if lists of target genes are in the filtered gene.expression matrix

intersect(ATG_markers$PhyCode, gene.ids) # Check

##  [1] "Cre09.g391245.v5.5" "Cre01.g045600.v5.5" "Cre02.g102350.v5.5"
##  [4] "Cre12.g510100.v5.5" "Cre14.g630907.v5.5" "Cre03.g165215.v5.5"
##  [7] "Cre16.g689650.v5.5" "Cre09.g391500.v5.5" "Cre12.g557000.v5.5"
## [10] "Cre16.g659000.v5.5" "Cre08.g377600.v5.5" "Cre10.g457550.v5.5"
## [13] "Cre07.g334700.v5.5" "Cre06.g290500.v5.5" "Cre01.g035500.v5.5"

intersect(Chloro_stress_markers$PhyCode, gene.ids) # Check

##  [1] "Cre13.g583550.v5.5" "Cre11.g468050.v5.5" "Cre06.g250100.v5.5"
##  [4] "Cre03.g145787.v5.5" "Cre14.g617450.v5.5" "Cre14.g617400.v5.5"
##  [7] "Cre02.g090850.v5.5" "Cre17.g696400.v5.5" "Cre12.g559150.v5.5"
## [10] "Cre07.g341600.v5.5" "Cre12.g507650.v5.5" "Cre01.g009900.v5.5"
## [13] "Cre10.g430400.v5.5" "Cre12.g520600.v5.5"

intersect(Chloro_redox_markers$PhyCode, gene.ids)  # Check

## [1] "Cre02.g087700.v5.5" "Cre06.g285150.v5.5" "Cre10.g458450.v5.5"
## [4] "Cre07.g325743.v5.5" "Cre16.g688550.v5.5" "Cre01.g054150.v5.5"
## [7] "Cre10.g422300.v5.5"

ATG Markers

bar.names_ATG <- ATG_markers$Name
bar.phycodes_ATG <- ATG_markers$PhyCode

wt_treatment_ATG <- as.data.frame(resWT) %>% rownames_to_column("GeneID")
wt_treatment_ATG <- wt_treatment_ATG[wt_treatment_ATG$GeneID %in% bar.phycodes_ATG,]
wt_treatment_ATG$GeneID <- ATG_markers$Name[match(wt_treatment_ATG$GeneID,ATG_markers$PhyCode)]
row.names(wt_treatment_ATG) <- wt_treatment_ATG$GeneID
wt_treatment_ATG_fold.change <- wt_treatment_ATG[-c(1:2,4:7)]
colnames(wt_treatment_ATG_fold.change) <- c("WT cerulenin")
wt_treatment_ATG_fold.change.matrix <- as.matrix(wt_treatment_ATG_fold.change)
wt_treatment_ATG$GeneID[wt_treatment_ATG$log2FoldChange >1 & wt_treatment_ATG$padj < 0.05]

## [1] "ATG7"   "ATG3"   "ATG14"  "ATG2"   "ATG101" "ATG8"

atg8_treatment_ATG <- as.data.frame(res_atg8) %>% rownames_to_column("GeneID")
atg8_treatment_ATG <- atg8_treatment_ATG[atg8_treatment_ATG$GeneID %in% bar.phycodes_ATG,]
atg8_treatment_ATG$GeneID <- ATG_markers$Name[match(atg8_treatment_ATG$GeneID,ATG_markers$PhyCode)]
row.names(atg8_treatment_ATG) <- atg8_treatment_ATG$GeneID
atg8_treatment_ATG_fold.change <- atg8_treatment_ATG[-c(1:2,4:7)]
colnames(atg8_treatment_ATG_fold.change) <- c("atg8 cerulenin")
atg8_treatment_ATG$GeneID[atg8_treatment_ATG$log2FoldChange >1 & atg8_treatment_ATG$padj < 0.05]

## [1] "VPS34"  "ATG7"   "ATG3"   "ATG14"  "ATG2"   "ATG101" "ATG8"   "ATG13"

atg8_treatment_ATG_fold.change.matrix <- as.matrix(atg8_treatment_ATG_fold.change)

wt_atg8_treatment_ATG_fold.change.matrix <- cbind(wt_treatment_ATG_fold.change.matrix,
                                              atg8_treatment_ATG_fold.change.matrix)
# Transpose matrix:

wt_atg8_treatment_ATG_fold.change.matrix <- t(wt_atg8_treatment_ATG_fold.change.matrix)

# Remove columns where for both genotypes target ATG genes have log2 fold change < 1
# and padj > 0.05

# Colums to remove: VPS15, ATG5, ATG12, ATG9, ATG5, ATG6, ATG1, ATG18, ATG4

wt_atg8_treatment_ATG_fold.change.matrix <- wt_atg8_treatment_ATG_fold.change.matrix[, -c(1,5,7:8,13:16)]
wt_atg8_treatment_ATG_fold.change.matrix<-wt_atg8_treatment_ATG_fold.change.matrix[,order(colnames(wt_atg8_treatment_ATG_fold.change.matrix), decreasing = TRUE)]

barplot(wt_atg8_treatment_ATG_fold.change.matrix, beside = TRUE,
        legend.text = TRUE, ylab = "log2 fold change", ylim = c(0,12),
        xpd = FALSE, cex.names = 0.5, col=colorRampPalette(rev(brewer.pal(2,"Blues")))(2))

## Warning in brewer.pal(2, "Blues"): minimal value for n is 3, returning requested palette with 3 different levels

Chloroplast chaperones and stress markers

bar.names_ChlStress <- Chloro_stress_markers$Name
bar.phycodes_ChlStress <- Chloro_stress_markers$PhyCode

wt_treatment_ChlStress <- as.data.frame(resWT) %>% rownames_to_column("GeneID")
wt_treatment_ChlStress <- wt_treatment_ChlStress[wt_treatment_ChlStress$GeneID %in% bar.phycodes_ChlStress,]
wt_treatment_ChlStress$GeneID <- Chloro_stress_markers$Name[match(wt_treatment_ChlStress$GeneID,
                                                                  Chloro_stress_markers$PhyCode)]
row.names(wt_treatment_ChlStress) <- wt_treatment_ChlStress$GeneID
wt_treatment_ChlStress_fold.change <- wt_treatment_ChlStress[-c(1:2,4:7)]
colnames(wt_treatment_ChlStress_fold.change) <- c("WT cerulenin")
wt_treatment_ChlStress_fold.change.matrix <- as.matrix(wt_treatment_ChlStress_fold.change)
wt_treatment_ChlStress$GeneID[wt_treatment_ChlStress$log2FoldChange >1 & wt_treatment_ChlStress$padj < 0.05]

## [1] "VIPP2"  "HSP22E" "VIPP1"  "HSP22F" "CLPB3"

atg8_treatment_ChlStress <- as.data.frame(res_atg8) %>% rownames_to_column("GeneID")
atg8_treatment_ChlStress <- atg8_treatment_ChlStress[atg8_treatment_ChlStress$GeneID %in% bar.phycodes_ChlStress,]
atg8_treatment_ChlStress$GeneID <- Chloro_stress_markers$Name[match(atg8_treatment_ChlStress$GeneID,
                                                                  Chloro_stress_markers$PhyCode)]
row.names(atg8_treatment_ChlStress) <- atg8_treatment_ChlStress$GeneID
atg8_treatment_ChlStress_fold.change <- atg8_treatment_ChlStress[-c(1:2,4:7)]
colnames(atg8_treatment_ChlStress_fold.change) <- c("atg8 cerulenin")
atg8_treatment_ChlStress$GeneID[atg8_treatment_ChlStress$log2FoldChange >1 & atg8_treatment_ChlStress$padj < 0.05]

## [1] "VIPP2"  "HSP22E" "VIPP1"  "HSP22F" "CLPB3"

atg8_treatment_ChlStress_fold.change.matrix <- as.matrix(atg8_treatment_ChlStress_fold.change)

wt_atg8_treatment_ChlStress_fold.change.matrix <- cbind(wt_treatment_ChlStress_fold.change.matrix,
                                              atg8_treatment_ChlStress_fold.change.matrix)
# Transpose matrix:

wt_atg8_treatment_ChlStress_fold.change.matrix <- t(wt_atg8_treatment_ChlStress_fold.change.matrix)

# Remove columns where for both genotypes target ATG genes have log2 fold change < 1
# and padj > 0.05

# Colums to remove: HSP22C, RPL37, CDJ2, CDJ1, CGE1, HSP70B, CDJ3, CLPS1, PRPS6, CLPS2

wt_atg8_treatment_ChlStress_fold.change.matrix <- wt_atg8_treatment_ChlStress_fold.change.matrix[, -c(2:3, 5:10,12,14)]
wt_atg8_treatment_ChlStress_fold.change.matrix<-wt_atg8_treatment_ChlStress_fold.change.matrix[,order(colnames(wt_atg8_treatment_ChlStress_fold.change.matrix), decreasing = TRUE)]

barplot(wt_atg8_treatment_ChlStress_fold.change.matrix, beside = TRUE,
        legend.text = TRUE, ylab = "log2 fold change", ylim = c(0,6),
        xpd = FALSE, cex.names = 0.5, col=colorRampPalette(rev(brewer.pal(2,"Blues")))(2))

## Warning in brewer.pal(2, "Blues"): minimal value for n is 3, returning requested palette with 3 different levels

Chloroplast redox markers

bar.names_ChlRed <- Chloro_redox_markers$Name
bar.phycodes_ChlRed <- Chloro_redox_markers$PhyCode

wt_treatment_ChlRed <- as.data.frame(resWT) %>% rownames_to_column("GeneID")
wt_treatment_ChlRed <- wt_treatment_ChlRed[wt_treatment_ChlRed$GeneID %in% bar.phycodes_ChlRed,]
wt_treatment_ChlRed$GeneID <- Chloro_redox_markers$Name[match(wt_treatment_ChlRed$GeneID,
                                                                  Chloro_redox_markers$PhyCode)]
row.names(wt_treatment_ChlRed) <- wt_treatment_ChlRed$GeneID
wt_treatment_ChlRed_fold.change <- wt_treatment_ChlRed[-c(1:2,4:7)]
colnames(wt_treatment_ChlRed_fold.change) <- c("WT cerulenin")
wt_treatment_ChlRed_fold.change.matrix <- as.matrix(wt_treatment_ChlRed_fold.change)
wt_treatment_ChlRed$GeneID[wt_treatment_ChlRed$log2FoldChange >1 & wt_treatment_ChlRed$padj < 0.05]

## [1] "GSTS1" "APX1"  "GPX5"

atg8_treatment_ChlRed <- as.data.frame(res_atg8) %>% rownames_to_column("GeneID")
atg8_treatment_ChlRed <- atg8_treatment_ChlRed[atg8_treatment_ChlRed$GeneID %in% bar.phycodes_ChlRed,]
atg8_treatment_ChlRed$GeneID <- Chloro_redox_markers$Name[match(atg8_treatment_ChlRed$GeneID,
                                                                  Chloro_redox_markers$PhyCode)]
row.names(atg8_treatment_ChlRed) <- atg8_treatment_ChlRed$GeneID
atg8_treatment_ChlRed_fold.change <- atg8_treatment_ChlRed[-c(1:2,4:7)]
colnames(atg8_treatment_ChlRed_fold.change) <- c("atg8 cerulenin")
atg8_treatment_ChlRed$GeneID[atg8_treatment_ChlRed$log2FoldChange >1 & atg8_treatment_ChlRed$padj < 0.05]

## [1] "GSTS1" "APX1"  "GPX5"

atg8_treatment_ChlRed_fold.change.matrix <- as.matrix(atg8_treatment_ChlRed_fold.change)

wt_atg8_treatment_ChlRed_fold.change.matrix <- cbind(wt_treatment_ChlRed_fold.change.matrix,
                                              atg8_treatment_ChlRed_fold.change.matrix)
# Transpose matrix:

wt_atg8_treatment_ChlRed_fold.change.matrix <- t(wt_atg8_treatment_ChlRed_fold.change.matrix)

# Remove columns where for both genotypes target ATG genes have log2 fold change < 1
# and padj > 0.05

# Colums to remove: APX2, GRX3, PRX6, NTRC1

wt_atg8_treatment_ChlRed_fold.change.matrix <- wt_atg8_treatment_ChlRed_fold.change.matrix[, -c(3:6)]
wt_atg8_treatment_ChlRed_fold.change.matrix<-wt_atg8_treatment_ChlRed_fold.change.matrix[,order(colnames(wt_atg8_treatment_ChlRed_fold.change.matrix), decreasing = TRUE)]


barplot(wt_atg8_treatment_ChlRed_fold.change.matrix, beside = TRUE,
        legend.text = TRUE, ylab = "log2 fold change", ylim = c(0,6),
        xpd = FALSE, cex.names = 0.5, col=colorRampPalette(rev(brewer.pal(2,"Blues")))(2))

## Warning in brewer.pal(2, "Blues"): minimal value for n is 3, returning requested palette with 3 different levels

Difference between wild-type and mutant 0h

Gene obtention

res_WT_atg8_0h = results(dds, contrast=c("genotype","atg8","wt"),
                         cooksCutoff = FALSE, independentFiltering = FALSE)
res_WT_atg8_0h

## log2 fold change (MLE): genotype atg8 vs wt 
## Wald test p-value: genotype atg8 vs wt 
## DataFrame with 13977 rows and 6 columns
##                      baseMean log2FoldChange     lfcSE      stat    pvalue
##                     <numeric>      <numeric> <numeric> <numeric> <numeric>
## Cre14.g619950.v5.5   307.8280     -1.0486329  0.719601 -1.457241  0.145050
## Cre01.g048900.v5.5    76.4131      0.2567535  1.830265  0.140282  0.888437
## Cre12.g543400.v5.5   428.9410      0.4253013  0.320785  1.325814  0.184901
## Cre16.g685837.v5.5   255.7288      0.0698002  0.456223  0.152996  0.878402
## Cre06.g278159.v5.5  9770.0557     -0.0796437  0.124401 -0.640217  0.522031
## ...                       ...            ...       ...       ...       ...
## Cre06.g278104.v5.5  889.54623      -0.604567  0.310855 -1.944851 0.0517930
## Cre17.g717150.v5.5  146.81853       0.994081  1.170782  0.849074 0.3958403
## Cre10.g443600.v5.5  170.10265       0.127679  1.079078  0.118322 0.9058126
## Cre14.g621400.v5.5    5.13639       6.786913  4.117852  1.648168 0.0993182
## Cre10.g441750.v5.5 1599.45474       0.269795  0.181866  1.483486 0.1379453
##                         padj
##                    <numeric>
## Cre14.g619950.v5.5  0.728842
## Cre01.g048900.v5.5  0.999313
## Cre12.g543400.v5.5  0.799370
## Cre16.g685837.v5.5  0.999313
## Cre06.g278159.v5.5  0.999313
## ...                      ...
## Cre06.g278104.v5.5  0.445209
## Cre17.g717150.v5.5  0.981999
## Cre10.g443600.v5.5  0.999313
## Cre14.g621400.v5.5  0.619202
## Cre10.g441750.v5.5  0.711462

log.fold.change_WT_atg8_0h <- res_WT_atg8_0h$log2FoldChange
q.value_WT_atg8_0h <- res_WT_atg8_0h$padj
names(log.fold.change_WT_atg8_0h) <- gene.ids
names(q.value_WT_atg8_0h) <- gene.ids
base.mean_WT_atg8_0h <- res_WT_atg8_0h$baseMean
names(base.mean_WT_atg8_0h) <- gene.ids

# Take a look at upregulated and downregulated genes: 

top_20_overexpressed_genes_WT_atg8_0h <- res_WT_atg8_0h[order(res_WT_atg8_0h$log2FoldChange, decreasing = T),]
kable(top_20_overexpressed_genes_WT_atg8_0h[1:20,-(3:4)])

	baseMean	log2FoldChange	pvalue	padj
Cre10.g466200.v5.5	15.235117	20.585945	0.0000013	0.0001480
Cre07.g331350.v5.5	19.178658	19.940809	0.0000005	0.0000721
Cre02.g143107.v5.5	7.114070	19.404147	0.0000051	0.0004478
Cre10.g465763.v5.5	17.445204	18.935435	0.0000086	0.0006709
Cre06.g250650.v5.5	13.336058	14.980359	0.0004177	0.0159081
Cre13.g567450.v5.5	47.816726	9.974608	0.0028980	0.0675091
Cre06.g310900.v5.5	12.697053	8.158493	0.0541851	0.4562322
Cre03.g167622.v5.5	13.099835	8.019118	0.0284711	0.3099224
Cre08.g358565.v5.5	7.390831	7.384139	0.0824190	0.5624718
Cre02.g120200.v5.5	7.613366	7.380495	0.0815050	0.5587032
Cre15.g635100.v5.5	16.735402	7.208932	0.0296900	0.3160524
Cre17.g713500.v5.5	5.698800	7.072735	0.0671212	0.5118263
Cre13.g569900.v5.5	12.543505	6.987735	0.1003772	0.6221609
Cre09.g404552.v5.5	7.466950	6.930540	0.1032099	0.6306568
Cre14.g621400.v5.5	5.136393	6.786913	0.0993182	0.6192016
Cre11.g480502.v5.5	4.849025	6.781913	0.0473071	0.4194994
Cre09.g395925.v5.5	10.571970	6.551579	0.1236302	0.6819178
Cre03.g144444.v5.5	8.399966	6.412211	0.1318948	0.6990198
Cre06.g265200.v5.5	3.559183	6.353395	0.1041611	0.6328170
Cre02.g091226.v5.5	16.082502	6.198403	0.1454064	0.7288417

top_20_overrepressed_genes_WT_atg8_0h <- res_WT_atg8_0h[order(res_WT_atg8_0h$log2FoldChange, decreasing = F),]
kable(top_20_overrepressed_genes_WT_atg8_0h[1:20,-(3:4)])

	baseMean	log2FoldChange	pvalue	padj
Cre07.g346900.v5.5	87.770943	-23.11461	0.0e+00	0.0000000
Cre08.g364650.v5.5	58.216483	-22.48245	0.0e+00	0.0000000
Cre02.g109900.v5.5	70.008000	-22.42246	0.0e+00	0.0000000
Cre03.g183150.v5.5	52.341644	-21.67106	0.0e+00	0.0000000
Cre10.g433650.v5.5	19.955260	-21.51465	1.0e-07	0.0000204
Cre02.g077951.v5.5	22.406577	-21.47873	1.0e-07	0.0000134
Cre03.g169150.v5.5	84.338786	-21.16178	0.0e+00	0.0000000
Cre02.g095141.v5.5	19.617552	-21.06563	3.0e-07	0.0000407
Cre13.g581700.v5.5	50.684675	-21.00907	0.0e+00	0.0000000
Cre10.g452500.v5.5	27.671211	-20.91627	5.0e-07	0.0000661
Cre06.g282826.v5.5	42.542451	-20.82930	0.0e+00	0.0000000
Cre12.g494600.v5.5	20.255885	-20.54698	3.0e-07	0.0000407
Cre07.g350050.v5.5	12.690641	-20.50715	1.4e-06	0.0001607
Cre01.g031700.v5.5	71.409961	-20.47323	0.0e+00	0.0000000
Cre03.g206900.v5.5	39.266900	-20.37765	0.0e+00	0.0000001
Cre09.g395732.v5.5	26.133934	-20.36151	1.0e-07	0.0000217
Cre06.g297049.v5.5	9.569606	-20.31073	1.2e-06	0.0001435
Cre17.g698950.v5.5	9.775073	-20.25779	1.9e-06	0.0002031
Cre03.g191650.v5.5	21.334812	-20.15767	0.0e+00	0.0000030
Cre16.g669350.v5.5	14.906514	-20.03924	9.0e-07	0.0001214

activated.genes.deseq2_WT_atg8_0h <- gene.ids[log.fold.change_WT_atg8_0h > 1 & q.value_WT_atg8_0h < 0.05]
activated.genes.deseq2_WT_atg8_0h <- activated.genes.deseq2_WT_atg8_0h[!is.na(activated.genes.deseq2_WT_atg8_0h)]
activated.genes.deseq2_WT_atg8_0h_good <- substr(activated.genes.deseq2_WT_atg8_0h,1,nchar(activated.genes.deseq2_WT_atg8_0h)-5)
write.table(activated.genes.deseq2_WT_atg8_0h_good,file="activated.genes.deseq2_WT_atg8_0h.txt",sep="\t",quote=F,row.names = F)

repressed.genes.deseq2_WT_atg8_0h <- gene.ids[log.fold.change_WT_atg8_0h < -1 & q.value_WT_atg8_0h < 0.05]
repressed.genes.deseq2_WT_atg8_0h <- repressed.genes.deseq2_WT_atg8_0h[!is.na(repressed.genes.deseq2_WT_atg8_0h)]
repressed.genes.deseq2_WT_atg8_0h_good <- substr(repressed.genes.deseq2_WT_atg8_0h,1,nchar(repressed.genes.deseq2_WT_atg8_0h)-5)
write.table(repressed.genes.deseq2_WT_atg8_0h_good,file="repressed.genes.deseq2_WT_atg8_0h.txt",sep="\t",quote=F,row.names = F)

length(activated.genes.deseq2_WT_atg8_0h_good)

## [1] 83

length(repressed.genes.deseq2_WT_atg8_0h_good)

## [1] 117

MA plot

plotMA(res_WT_atg8_0h,alpha=0.05, ylim=c(-30,30),colSig="grey")
abline(h=c(-1,1), col="black", lwd=1)
points(x = base.mean_WT_atg8_0h[activated.genes.deseq2_WT_atg8_0h],
       y = log.fold.change_WT_atg8_0h[activated.genes.deseq2_WT_atg8_0h],col="red",cex=0.2,pch=19)
points(x = base.mean_WT_atg8_0h[repressed.genes.deseq2_WT_atg8_0h],
       y = log.fold.change_WT_atg8_0h[repressed.genes.deseq2_WT_atg8_0h],col="blue",cex=0.2,pch=19)

text(x =100000, y = 15, label = length(activated.genes.deseq2_WT_atg8_0h_good),
     col = "red",   # Color del texto
     font = 2,      # Estilo
     cex = 1.5)     # Tamaño
text(x =100000, y = -15, label = length(repressed.genes.deseq2_WT_atg8_0h_good),
     col = "blue",   # Color del texto
     font = 2,      # Estilo
     cex = 1.5)     # Tamaño

legend("bottomright", legend = c("Overexpressed", "Repressed"),
       lwd = 3, col = c("red", "blue"))

Difference between wild-type and mutant 4h

Gene obtention

res_WT_atg8_4h = results(dds, list( c("genotype_atg8_vs_wt","genotypeatg8.condition4h")),
                                    cooksCutoff = FALSE, independentFiltering = FALSE )
res_WT_atg8_4h

## log2 fold change (MLE): genotype_atg8_vs_wt+genotypeatg8.condition4h effect 
## Wald test p-value: genotype_atg8_vs_wt+genotypeatg8.condition4h effect 
## DataFrame with 13977 rows and 6 columns
##                      baseMean log2FoldChange     lfcSE      stat     pvalue
##                     <numeric>      <numeric> <numeric> <numeric>  <numeric>
## Cre14.g619950.v5.5   307.8280      -1.094574  0.719128 -1.522086 0.12798748
## Cre01.g048900.v5.5    76.4131      -0.476455  1.829920 -0.260369 0.79457914
## Cre12.g543400.v5.5   428.9410       0.631098  0.320734  1.967667 0.04910641
## Cre16.g685837.v5.5   255.7288       0.130937  0.454228  0.288264 0.77314498
## Cre06.g278159.v5.5  9770.0557       0.361482  0.124145  2.911778 0.00359378
## ...                       ...            ...       ...       ...        ...
## Cre06.g278104.v5.5  889.54623       0.228314  0.310320  0.735738  0.4618902
## Cre17.g717150.v5.5  146.81853      -0.644081  1.173279 -0.548958  0.5830341
## Cre10.g443600.v5.5  170.10265       1.887998  1.082569  1.743998  0.0811594
## Cre14.g621400.v5.5    5.13639       4.605473  4.141349  1.112071  0.2661077
## Cre10.g441750.v5.5 1599.45474       0.277861  0.181312  1.532499  0.1253994
##                         padj
##                    <numeric>
## Cre14.g619950.v5.5 0.5198724
## Cre01.g048900.v5.5 0.9823882
## Cre12.g543400.v5.5 0.3269940
## Cre16.g685837.v5.5 0.9793068
## Cre06.g278159.v5.5 0.0577359
## ...                      ...
## Cre06.g278104.v5.5  0.879075
## Cre17.g717150.v5.5  0.930899
## Cre10.g443600.v5.5  0.420310
## Cre14.g621400.v5.5  0.733029
## Cre10.g441750.v5.5  0.515712

log.fold.change_WT_atg8_4h <- res_WT_atg8_4h$log2FoldChange
q.value_WT_atg8_4h <- res_WT_atg8_4h$padj
names(log.fold.change_WT_atg8_4h) <- gene.ids
names(q.value_WT_atg8_4h) <- gene.ids
base.mean_WT_atg8_4h <- res_WT_atg8_4h$baseMean
names(base.mean_WT_atg8_4h) <- gene.ids

# Take a look at upregulated and downregulated genes: 

top_20_overexpressed_genes_WT_atg8_4h <- res_WT_atg8_4h[order(res_WT_atg8_4h$log2FoldChange, decreasing = T),]
kable(top_20_overexpressed_genes_WT_atg8_4h[1:20,-(3:4)])

	baseMean	log2FoldChange	pvalue	padj
Cre17.g712700.v5.5	42.409362	23.20189	0e+00	0.00e+00
Cre16.g658400.v5.5	59.759725	23.01070	0e+00	0.00e+00
Cre12.g524450.v5.5	32.024953	21.89046	0e+00	4.00e-07
Cre13.g588250.v5.5	49.571334	21.80803	0e+00	0.00e+00
Cre16.g687966.v5.5	33.784426	21.79817	0e+00	0.00e+00
Cre03.g143827.v5.5	66.501562	21.66314	0e+00	0.00e+00
Cre15.g636600.v5.5	58.721594	21.62000	0e+00	0.00e+00
Cre07.g336500.v5.5	36.496724	21.41570	0e+00	4.00e-07
Cre01.g043250.v5.5	102.227854	21.36041	0e+00	0.00e+00
Cre03.g195410.v5.5	16.175576	21.30459	0e+00	5.00e-07
Cre12.g524350.v5.5	18.072167	21.27414	6e-07	3.66e-05
Cre03.g179400.v5.5	39.299172	21.24146	0e+00	0.00e+00
Cre11.g467532.v5.5	33.222102	21.23809	0e+00	0.00e+00
Cre05.g240100.v5.5	25.335507	21.21101	0e+00	3.60e-06
Cre05.g239000.v5.5	46.353687	21.15574	0e+00	0.00e+00
Cre06.g292600.v5.5	15.273678	21.14644	7e-07	4.17e-05
Cre01.g034380.v5.5	38.571871	21.10793	0e+00	0.00e+00
Cre10.g422350.v5.5	8.728187	21.10641	7e-07	4.32e-05
Cre09.g402050.v5.5	54.881325	21.10392	0e+00	0.00e+00
Cre14.g629840.v5.5	37.019186	21.03051	0e+00	0.00e+00

top_20_overrepressed_genes_WT_atg8_4h <- res_WT_atg8_4h[order(res_WT_atg8_4h$log2FoldChange, decreasing = F),]
kable(top_20_overrepressed_genes_WT_atg8_4h[1:20,-(3:4)])

	baseMean	log2FoldChange	pvalue	padj
Cre02.g077951.v5.5	22.40658	-27.39922	0e+00	0.00e+00
Cre10.g433650.v5.5	19.95526	-26.86531	0e+00	0.00e+00
Cre07.g350050.v5.5	12.69064	-26.22043	0e+00	1.00e-07
Cre02.g095141.v5.5	19.61755	-26.18287	0e+00	0.00e+00
Cre07.g315050.v5.5	88.92973	-23.62099	0e+00	0.00e+00
Cre16.g684861.v5.5	36.98465	-22.61599	0e+00	0.00e+00
Cre07.g327079.v5.5	23.64919	-22.52874	1e-07	6.40e-06
Cre13.g592350.v5.5	49.26531	-22.51609	0e+00	0.00e+00
Cre14.g623300.v5.5	27.93290	-22.39287	0e+00	0.00e+00
Cre13.g604550.v5.5	56.05740	-22.25126	0e+00	0.00e+00
Cre02.g118450.v5.5	30.20969	-22.24960	0e+00	4.50e-06
Cre03.g205800.v5.5	69.35132	-22.19388	2e-07	1.33e-05
Cre12.g509600.v5.5	38.66322	-22.15181	0e+00	0.00e+00
Cre06.g278168.v5.5	57.62603	-22.07722	0e+00	0.00e+00
Cre03.g153250.v5.5	28.01849	-22.07050	1e-07	9.40e-06
Cre07.g314700.v5.5	37.60714	-21.95988	0e+00	0.00e+00
Cre16.g678997.v5.5	30.44633	-21.93106	0e+00	0.00e+00
Cre06.g294250.v5.5	123.06077	-21.90895	0e+00	0.00e+00
Cre15.g636250.v5.5	24.67298	-21.89186	0e+00	1.80e-06
Cre01.g010832.v5.5	37.26769	-21.87643	0e+00	0.00e+00

activated.genes.deseq2_WT_atg8_4h <- gene.ids[log.fold.change_WT_atg8_4h > 1 & q.value_WT_atg8_4h < 0.05]
activated.genes.deseq2_WT_atg8_4h <- activated.genes.deseq2_WT_atg8_4h[!is.na(activated.genes.deseq2_WT_atg8_4h)]
activated.genes.deseq2_WT_atg8_4h_good <- substr(activated.genes.deseq2_WT_atg8_4h,1,nchar(activated.genes.deseq2_WT_atg8_4h)-5)
write.table(activated.genes.deseq2_WT_atg8_4h_good,file="activated.genes.deseq2_WT_atg8_4h.txt",sep="\t",quote=F,row.names = F)


repressed.genes.deseq2_WT_atg8_4h <- gene.ids[log.fold.change_WT_atg8_4h < -1 & q.value_WT_atg8_4h < 0.05]
repressed.genes.deseq2_WT_atg8_4h <- repressed.genes.deseq2_WT_atg8_4h[!is.na(repressed.genes.deseq2_WT_atg8_4h)]
repressed.genes.deseq2_WT_atg8_4h_good <- substr(repressed.genes.deseq2_WT_atg8_4h,1,nchar(repressed.genes.deseq2_WT_atg8_4h)-5)
write.table(repressed.genes.deseq2_WT_atg8_4h_good,file="repressed.genes.deseq2_WT_atg8_4h.txt",sep="\t",quote=F,row.names = F)


length(activated.genes.deseq2_WT_atg8_4h)

## [1] 193

length(repressed.genes.deseq2_WT_atg8_4h)

## [1] 187

MA plot

plotMA(res_WT_atg8_4h, ylim=c(-30,30),colSig="grey")
abline(h=c(-1,1), col="black", lwd=1)
points(x = base.mean_WT_atg8_4h[activated.genes.deseq2_WT_atg8_4h],
       y = log.fold.change_WT_atg8_4h[activated.genes.deseq2_WT_atg8_4h],col="red",cex=0.2,pch=19)
points(x = base.mean_WT_atg8_4h[repressed.genes.deseq2_WT_atg8_4h],
       y = log.fold.change_WT_atg8_4h[repressed.genes.deseq2_WT_atg8_4h],col="blue",cex=0.2,pch=19)

text(x =100000, y = 15, label = length(activated.genes.deseq2_WT_atg8_4h_good),
     col = "red",   # Color del texto
     font = 2,      # Estilo
     cex = 1.5)     # Tamaño
text(x =100000, y = -15, label = length(repressed.genes.deseq2_WT_atg8_4h_good),
     col = "blue",   # Color del texto
     font = 2,      # Estilo
     cex = 1.5)     # Tamaño

legend("bottomright", legend = c("Overexpressed", "Repressed"),
       lwd = 3, col = c("red", "blue"))

Different response for both genotypes (the interaction term)

Gene obtention

res_interaction = results(dds, name="genotypeatg8.condition4h",
                          cooksCutoff = FALSE, independentFiltering = FALSE)
res_interaction

## log2 fold change (MLE): genotypeatg8.condition4h 
## Wald test p-value: genotypeatg8.condition4h 
## DataFrame with 13977 rows and 6 columns
##                      baseMean log2FoldChange     lfcSE       stat    pvalue
##                     <numeric>      <numeric> <numeric>  <numeric> <numeric>
## Cre14.g619950.v5.5   307.8280     -0.0459415  1.017335 -0.0451587 0.9639808
## Cre01.g048900.v5.5    76.4131     -0.7332081  2.588131 -0.2832963 0.7769497
## Cre12.g543400.v5.5   428.9410      0.2057964  0.453622  0.4536733 0.6500640
## Cre16.g685837.v5.5   255.7288      0.0611372  0.643788  0.0949648 0.9243428
## Cre06.g278159.v5.5  9770.0557      0.4411260  0.175749  2.5099837 0.0120737
## ...                       ...            ...       ...        ...       ...
## Cre06.g278104.v5.5  889.54623     0.83288039  0.439237   1.896199 0.0579338
## Cre17.g717150.v5.5  146.81853    -1.63816151  1.657499  -0.988333 0.3229896
## Cre10.g443600.v5.5  170.10265     1.76031972  1.528515   1.151653 0.2494635
## Cre14.g621400.v5.5    5.13639    -2.18143969  5.840032  -0.373532 0.7087524
## Cre10.g441750.v5.5 1599.45474     0.00806577  0.256806   0.031408 0.9749441
##                         padj
##                    <numeric>
## Cre14.g619950.v5.5  0.999419
## Cre01.g048900.v5.5  0.999419
## Cre12.g543400.v5.5  0.999419
## Cre16.g685837.v5.5  0.999419
## Cre06.g278159.v5.5  0.366382
## ...                      ...
## Cre06.g278104.v5.5  0.772653
## Cre17.g717150.v5.5  0.999419
## Cre10.g443600.v5.5  0.999419
## Cre14.g621400.v5.5  0.999419
## Cre10.g441750.v5.5  0.999419

log.fold.change_interaction <- res_interaction$log2FoldChange
q.value_interaction <- res_interaction$padj
names(log.fold.change_interaction) <- gene.ids
names(q.value_interaction) <- gene.ids
base.mean_interaction <- res_interaction$baseMean
names(base.mean_interaction) <- gene.ids

# Take a look at upregulated and downregulated genes: 

top_20_overexpressed_genes_interaction <- res_interaction[order(res_interaction$log2FoldChange, decreasing = T),]
kable(top_20_overexpressed_genes_interaction[1:20,-(3:4)])

	baseMean	log2FoldChange	pvalue	padj
Cre04.g214769.v5.5	8.578910	30.00000	0.0000006	0.0001753
Cre07.g331114.v5.5	14.030683	30.00000	0.0000006	0.0001753
Cre03.g196000.v5.5	10.329530	30.00000	0.0000006	0.0001753
Cre04.g213950.v5.5	7.990901	30.00000	0.0000006	0.0001753
Cre06.g292450.v5.5	8.428540	30.00000	0.0000006	0.0001753
Cre11.g481104.v5.5	17.175085	30.00000	0.0000003	0.0001134
Cre04.g219750.v5.5	10.123388	30.00000	0.0000006	0.0001753
Cre02.g087500.v5.5	10.921086	30.00000	0.0000006	0.0001753
Cre14.g633900.v5.5	24.339891	30.00000	0.0000001	0.0000342
Cre02.g099100.v5.5	33.227640	30.00000	0.0000000	0.0000036
Cre16.g658300.v5.5	46.217102	23.52922	0.0000000	0.0000006
Cre03.g183150.v5.5	52.341644	22.83219	0.0000000	0.0000245
Cre03.g206900.v5.5	39.266900	22.28228	0.0000005	0.0001753
Cre07.g319600.v5.5	9.767238	22.09547	0.0001150	0.0143253
Cre02.g098550.v5.5	10.725460	22.02765	0.0001957	0.0199689
Cre17.g699550.v5.5	11.544327	22.01875	0.0001931	0.0199689
Cre07.g346900.v5.5	87.770943	21.99156	0.0000000	0.0000000
Cre07.g330150.v5.5	6.660717	21.94332	0.0002175	0.0218705
Cre03.g179400.v5.5	39.299172	21.70125	0.0000006	0.0001753
Cre16.g658400.v5.5	59.759725	21.68779	0.0000001	0.0000285

top_20_overrepressed_genes_interaction <- res_interaction[order(res_interaction$log2FoldChange, decreasing = F),]
kable(top_20_overrepressed_genes_interaction[1:20,-(3:4)])

	baseMean	log2FoldChange	pvalue	padj
Cre06.g250650.v5.5	13.336058	-30.00000	0.0000006	0.0001753
Cre16.g691351.v5.5	18.498013	-24.96833	0.0000052	0.0012079
Cre07.g315050.v5.5	88.929727	-24.78043	0.0000000	0.0000000
Cre02.g118450.v5.5	30.209685	-24.24031	0.0000181	0.0034299
Cre13.g592350.v5.5	49.265306	-24.11483	0.0000000	0.0000053
Cre13.g604550.v5.5	56.057399	-23.90740	0.0000000	0.0000001
Cre16.g678997.v5.5	30.446334	-23.88634	0.0000003	0.0001089
Cre16.g695200.v5.5	28.693675	-23.54711	0.0000408	0.0066344
Cre12.g544327.v5.5	9.507146	-23.34612	0.0000802	0.0112272
Cre09.g413000.v5.5	10.374182	-23.30969	0.0000816	0.0112272
Cre03.g153250.v5.5	28.018485	-23.29247	0.0000548	0.0086987
Cre16.g684861.v5.5	36.984651	-23.08733	0.0000000	0.0000077
Cre12.g486450.v5.5	11.227031	-22.92555	0.0001045	0.0133968
Cre08.g382825.v5.5	12.828425	-22.85710	0.0001092	0.0138783
Cre14.g609150.v5.5	11.132919	-22.63246	0.0001300	0.0152511
Cre03.g151300.v5.5	13.711523	-22.62959	0.0001261	0.0150679
Cre03.g205800.v5.5	69.351320	-22.56410	0.0001282	0.0151822
Cre10.g466200.v5.5	15.235117	-22.51100	0.0001425	0.0161548
Cre17.g736350.v5.5	15.843663	-22.49394	0.0001378	0.0159119
Cre04.g215001.v5.5	11.408657	-22.48296	0.0001406	0.0161085

activated.genes.deseq2_interaction <- gene.ids[log.fold.change_interaction > 1 & q.value_interaction < 0.05]
activated.genes.deseq2_interaction <- activated.genes.deseq2_interaction[!is.na(activated.genes.deseq2_interaction)]
activated.genes.deseq2_interaction_good <- substr(activated.genes.deseq2_interaction,1,nchar(activated.genes.deseq2_interaction)-5)
write.table(activated.genes.deseq2_interaction_good,file="activated.genes.deseq2_interaction.txt",sep="\t",quote=F,row.names = F)


repressed.genes.deseq2_interaction <- gene.ids[log.fold.change_interaction < -1 & q.value_interaction < 0.05]
repressed.genes.deseq2_interaction <- repressed.genes.deseq2_interaction[!is.na(repressed.genes.deseq2_interaction)]
repressed.genes.deseq2_interaction_good <- substr(repressed.genes.deseq2_interaction,1,nchar(repressed.genes.deseq2_interaction)-5)
write.table(repressed.genes.deseq2_interaction_good,file="repressed.genes.deseq2_interaction.txt",sep="\t",quote=F,row.names = F)


length(activated.genes.deseq2_interaction)

## [1] 103

length(repressed.genes.deseq2_interaction)

## [1] 87

MA plot

plotMA(res_interaction, ylim=c(-30,30),colSig="grey")
abline(h=c(-1,1), col="black", lwd=1)
points(x = base.mean_interaction[activated.genes.deseq2_interaction],
       y = log.fold.change_interaction[activated.genes.deseq2_interaction],col="red",cex=0.2,pch=19)
points(x = base.mean_interaction[repressed.genes.deseq2_interaction],
       y = log.fold.change_interaction[repressed.genes.deseq2_interaction],col="blue",cex=0.2,pch=19)

text(x =100000, y = 15, label = length(activated.genes.deseq2_interaction_good),
     col = "red",   # Color del texto
     font = 2,      # Estilo
     cex = 1.5)     # Tamaño
text(x =100000, y = -15, label = length(repressed.genes.deseq2_interaction_good),
     col = "blue",   # Color del texto
     font = 2,      # Estilo
     cex = 1.5)     # Tamaño

legend("bottomright", legend = c("Overexpressed", "Repressed"),
       lwd = 3, col = c("red", "blue"))

Heatmap

interaction_sigs <- res_interaction[res_interaction$padj < 0.05,]
interaction_df <- as.data.frame(interaction_sigs)
interaction_dftop <- interaction_df[interaction_df$baseMean > 50 & 
                                    abs(interaction_df$log2FoldChange) > 1,]
interaction_dftop <- interaction_dftop[order(interaction_dftop$log2FoldChange, decreasing=TRUE),]
mat_int <- assay(vsd)[rownames(interaction_dftop), rownames(colData(vsd))]
colnames(mat_int) <- rownames(colData(vsd))
base_mean_int <- rowMeans(mat_int)
mat_int_scaled <- t(apply(mat_int, 1,scale))
colnames(mat_int_scaled) <- colnames(mat_int) # Original code, changed sample names in final figure
num_keep <- 25
rows_keep <- c(seq(1:num_keep), seq((nrow(mat_int_scaled)-num_keep),nrow(mat_int_scaled)))
l2_val <- as.matrix(interaction_dftop[rows_keep,]$log2FoldChange)
colnames(l2_val) <- "logFC"
mean <- as.matrix(interaction_dftop[rows_keep,]$baseMean)
colnames(mean) <- "AveExp"

# maps values between b/w/r for min and max l2 values
col_logFC <- colorRamp2(c(min(l2_val),0, max(l2_val)),c("blue", "white","red"))

# maps between 0% quantile, and 75% quantile of mean values
col_AveExpr <- colorRamp2(c(quantile(mean)[1],quantile(mean)[4]),c("white", "red"))

ha <- HeatmapAnnotation(summary = anno_summary(gp=gpar(fill=2), height =unit(2, "cm")))

h1 <- Heatmap(mat_int_scaled[rows_keep,],cluster_rows = F, column_labels = colnames(mat_int_scaled), name="Z-score", cluster_columns = T)

h2 <- Heatmap(l2_val, row_labels = rownames(interaction_dftop[rows_keep,]), cluster_rows = F, name= "logFC", top_annotation = ha, col= col_logFC, cell_fun = function(j,i,x,y,w,h,col){
  grid.text(round(l2_val[i,j],2),x,y)
})

h3 <- Heatmap(mean, row_labels = rownames(interaction_dftop[rows_keep,]), cluster_rows = F, name= "AveExpr", col = col_AveExpr, cell_fun = function(j,i,x,y,w,h,col){
  grid.text(round(mean[i,j],2),x,y)
})

h <- h1+h2+h3
h

Are there any common genes between WT and atg8 the activated and/or repressed genes with treatment?

Venn Diagrams

ven <- venndetail(list(WT_cerulenin= activated.genes.deseq2_WT_good, atg8_cerulenin=activated.genes.deseq2_atg8_good))
plot(ven, order=FALSE)

dplot(ven, order = TRUE, textsize = 2)

plot(ven, type = "upset")

ven2 <- venndetail(list(WT_cerulenin= repressed.genes.deseq2_WT_good, atg8_cerulenin= repressed.genes.deseq2_atg8_good))
plot(ven2, cat.cex=1, order=FALSE)

dplot(ven2, order = TRUE, textsize = 2)

plot(ven2, type = "upset")

Lists of shared genes

shared.activated <- result(ven)
shared.activated_cerulenin <- shared.activated %>% filter(shared.activated$Subset=="Shared")
shared.repressed <- result(ven2)
shared.repressed_cerulenin <- shared.repressed %>% filter(shared.repressed$Subset=="Shared")

Are there any common genes between WT and atg8 activated and/or repressed genes with and without treatment?

Venn Diagrams

ven3 <- venndetail(list(WTvsatg8_0h= activated.genes.deseq2_WT_atg8_0h_good, WTvsatg8_4h=activated.genes.deseq2_WT_atg8_4h_good))
plot(ven3)

dplot(ven3, order = TRUE, textsize = 2)

plot(ven3, type = "upset")

ven4 <- venndetail(list(WTvsatg8_0h= repressed.genes.deseq2_WT_atg8_0h_good, WTvsatg8_4h= repressed.genes.deseq2_WT_atg8_4h_good))
plot(ven4)

dplot(ven4, order = TRUE, textsize = 2)

plot(ven4, type = "upset")

Lists of shared genes

shared.activated2 <- result(ven3)
shared.activated_0h <- shared.activated2 %>% filter(shared.activated2$Subset=="Shared")
shared.repressed2 <- result(ven4)
shared.repressed_4h <- shared.repressed2 %>% filter(shared.repressed2$Subset=="Shared")